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1. Introduction 

A longstanding issue in general relativity has been to find the general behavior of the 
singularities that form in gravitational collapse. As far back as 1971, a conjecture on 
this subject was formulated by Belinskii, Khalatnikov and Lifshitz (BKL).[l] However, 
BKL provided no means of verifying their conjecture. To remedy this situation, Berger 
and Moncrief and their collaborators[2l[3],|l],[5j[6] i n the 1990s embarked on a program of 
examining the nature of singularities by numerical simulations: choose initial data whose 
evolution will become singular and evolve it numerically to see whether the singularity 
that forms is of the type conjectured by BKL. 

As is often the case in numerical relativity, the Berger-Moncrief plan was to first 
treat spacetimes with symmetry and then eventually proceed to the general case of 
spacetimes with no symmetry. Cases with symmetry are often easier to treat numerically 
and require a great deal less in the way of computational resources (or equivalently, with 
a fixed amount of computational resources assuming symmetry allows one to obtain 
much more resolution). However, spacetimes with a particular symmetry may exhibit 
behavior different from that of the general spacetime, so eventually one needs to tackle 
the general problem. 

Many of the issues involved in numerical simulations of singularities appear in the 
first and simplest models studied by Berger and Moncrief: the Gowdy spacetimes. [7] 
These are vacuum spacetimes with two spacelike, commuting Killing fields that are 
orthogonally transitive. Other systems with symmetry treated in this research program 
include the general vacuum T 2 symmetric spacetime (no orthogonal transitivity) and 
the vacuum U(l) symmetric spacetimes (only a single Killing field). Nonetheless, in all 
these cases the method made use in an essential way of the symmetry, and there was no 
clear way to generalize these methods to the case of no symmetry: a different method 
was needed. 

This different method was provided by Uggla et al[8], using a set of scale invariant 
variables in a tetrad formalism. Though this sort of method had originally been 
developed to treat homogeneous spacetimes, it was generalized without too much 
difficulty to the case with no symmetry, and it turns out to be suitable for a numerical 
treatment. 

In this paper, I will tell the begining and end of the story while omitting the middle. 
Section 2 will cover the Gowdy spacetimes and what has been learned about them 
from numerical simulations. Section 3 treats the Uggla et al system and its numerical 
implementation. Section 4 will consider open problems and areas for improvement in 
the numerical simulation of singularities. 

2. Gowdy Spacetimes 

The Gowdy spacetimes on T 3 x R take the form 

ds 2 = e ix+T)/2 (-e- 2T dr 2 + dx 2 ) + e- T [e p (dy + Qdzf + e~ p dz 2 } (1) 
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where the metric functions P, Q and A depend only on the time coordinate r and the 
spatial coordinate x. Thus (d/dy) a and (d/dz) a are Killing fields. The T 3 spatial 
topology is imposed by having < x,y,z < 2tt and having P, Q and A be periodic 
functions of x. 

The time coordinate r has an invariant geometric meaning, since the area of the 
orbits of the symmetry group is 47r 2 e~ T . This area approaches zero as the singularity 
is approached, and thus in terms of coordinates the singularity is at r = oo (though it 
can certainly be reached by observers in finite proper time). This property of the time 
coordinate, that the singularity is reached at infinite time, is very helpful for numerical 
simulations. A numerical simulation must stop at the earliest time when any part of 
a time slice reaches the singularity. Thus an ideal time coordinate is one in which no 
constant time slice ever hits the singularity, but in which the slices of constant time 
approach arbitrarily close to the singularity as the time goes to infinity. In this way, 
the behavior of the spacetime as the singularity is approached can be read off as the 
limiting behavior of the spacetime at large values of the time coordinate. 

For spacetimes of the form in eqn. ([1]) the vacuum Einstein equations become 

P,rr = e 2P Q 2 T + e~^P. xx - e^Q% (2) 

Q,tt = - 2P T Q tT + e~ 2T (Q :XX + 2P >X Q >X ) (3) 

as well as equations (essentially the usual constraint equations) that determine A once 
P and Q are known. Here the notation for derivatives is that a = 8/ da. 

One part of the BKL conjecture is that as the singularity is approached, terms in 
the Einstein field equations containing spatial derivatives become negligible compared 
to those that contain time derivatives. A glance at equations (TSJ) and ([3]) indicates how 
that might come about: in these equations, all the terms that involve spatial derivatives 
contain a factor of e~ 2r . As r — > oo one might therefore expect these terms to be 
negligible. Neglecting those terms allows one to solve equations (jSJ) and (J3j) in closed 
form and leads to the expected asymptotic behavior for large r 

P — > p(x) + tv(x) (4) 

Q - q(x) (5) 

for some functions]? (a:), v(x) and q(x). However, a second glance at equation <^ reveals 
that things are not quite so simple: the last term in that equation contains not the 
exponential of — 2r but instead the exponential of 2(P — r). Thus if P grows faster than 
r this term might not be negligible. One can resolve this question by doing a numerical 
simulation as was done in[2]. There are many ways to perform such a simulation, and 
in fact the method of [2] uses a sophisticated technique called symplectic integration. 
However, one can just as easily use the more standard numerical relativity technique 
of using centered differences for spatial derivatives, and the Iterative Crank-Nicholson 
(ICN) method for time derivatives (after first rewriting the system as one that is first 
order in time). The results of such a simulation are given in figures [T] and [2j 

Here the initial data at r = are P = 0, P T = 5cosx, Q = cosx, Q T = and 
the system has been evolved to r = 10. Note the presence of narrow features in these 
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Figure 1. Plot of P vs x for vq = 5 at r = 10. Here < x < ir 
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graphs that have been termed "spikes." Away from the spikes, the asymptotic form of 
the solution is that of equations (TjJ and (jSJ) with v(x) < 1. However, at the upward 
pointing spikes in P, we have that P grows faster than r and that the spike becomes 
ever narrower as the singularity is approached. Thus the simulations do show that as the 
singularity is approached the dynamics becomes one where the spatial derivatives can 
be neglected; however, they also show that there are special points where the dynamics 
is different from the neighboring points. This is a challenge for the numerics even for 
a 1+1 dimensional simulation. One can use adaptive mesh refinement to resolve the 
spikes [9 J or simply use a very large number of spatial points; but the spikes narrow 
exponentially in time and thus will eventually become too narrow for the numerical 
method to resolve them. 

3. scale invariant tetrad systems 

The Gowdy spacetimes have two features that are very convenient for the study of 
singularities: the first is the time coordinate that goes to infinity as the singularity 
is approached. The second is that even though the metric components themselves 
are exponential functions of the time coordinate, the quantities P and Q that are 
evolved by the numerics have much more mild behavior, growing no more than linearly 
with time. One would like to have both of these features for a system that has no 
symmetries whatsoever. It turns out that a convenient approach came from the study 
of homogeneous spacetimes using the tetrad formalism. In a homogeneous spacetime, 
each time slice has constant mean curvature, and while various geometric quantities 
blow up as the singularity is approached, they do so as particular powers of the mean 
curvature. Thus, one can introduce scale invariant variables by dividing each geometric 
quantity in the equations of motion by the appropriate power of mean curvature and 
then rewriting the field equations in terms of these scale invariant quantities. 

Not only does this work for homogeneous spacetimes; but it can also be generalized 
to the case of no symmetry at all, as was done by Uggla et al[8]. The system of reference 
[8] is the following: the spacetime is described in terms of a coordinate system (t, x l ) 
and a tetrad (e , e a ) where both the spatial coordinate index i and the spatial tetrad 
index a go from 1 to 3. Choose e to be hypersurface orthogonal with the relation 
between tetrad and coordinates of the form eo = N~ x dt and e a = e a l di where N is the 
lapse and the shift is chosen to be zero. Choose the spatial frame {e Q } to be Fermi 
propagated along the integral curves of eo- The commutators of the tetrad components 
are decomposed as follows: 



where n a ^ is symmetric, and a"? is symmetric and trace free. Scale invariant variables 
are defined as follows: {d ,d a } = {e ,e a }/H 




(6) 
(7) 



{E a \ Sq/3, A a , N a/3 } = {e a \a a (3 } a a } n al3 } / H 



(8) 
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q + 1 = — d In H and r a = —d a In H. 

Finally choose the lapse to be N = H~ l . In the conventions of [8] the singularity is 
in the negative t direction and one expects it to be approached as t — > — oo. The relation 
between scale invariant frame derivatives and coordinate derivatives is d = d t and 
d a = E a l di. From the vacuum Einstein equations one obtains the following evolution 
equations: 



d t E a l 
d t r a 
d t A a 
<9 t £ Q/3 

d t N ap 
d t q 



Fjrp + d a q 
F a l3 A 13 + \dpY?P 

(q - 2)S Q/3 - 2N <a 7 N p> ^ + NJN <af5> + d <a r p> 
q«* A P> + 2r <a A fi> + e lS<a {d 1 - 2A J )N? > S 

2(q - 2) + | (2A a - r a ) d a - \d a 
\d a r a + \A a r a + lr p d a E a P - 2Z af3 W ( 



q 

a/3 



(9) 
(10) 

(11) 

(12) 
(13) 

(14) 



Here angle brackets denote the symmetric trace- free part, and F a p and W a p are given 
by 



F a p 



q5 a p — S Q/ g 

§7V Q7 iV - \N\N aP + \d a A(. 



d a r p - \e^ a (d 1 - 2A,) N f 
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(15) 



(16) 



17) 



In addition to the evolution equations, the variables are subject to the vanishing of the 
following constraints: 

{Ccomfap = 2(d[ a — T[ a — A[ a )Ep] 1 — t a psN 5 "' ' Ej 

C G = 1 + \{2d a - 2r a - 3A a )A a - \N a pN^ 

= dpE a/3 + 2r a - S^r^ - 3ApE a/3 - e a ^N 08 ^ S 



n — iV^V „ _|_ I/) r a — 1 A r a 



{Cm 



(18) 
(19) 
(20) 
(21) 
(22) 



= (dp - rp)(N af3 + e a ^A y ) - 2ApN a ? 
= [e a ^(dp - Ap) - AHr 7 

In fact, this set of equations yields a class of evolution equations, since one can always 
add any of the constraint equations to the right hand side of any of the evolution 
equations. In fact, the system as written is not well posed, but it becomes well posed if 
one adds d times equation (fl9l to the right hand side of equation (TTTj) where d is any 
number less then — 1/2.[T0] 

The properties of the time coordinate come about in the choice of lapse: iV = H^ 1 
where 3H is the mean curvature. This is sometimes known as inverse mean curvature 
flow: a system that has been studied in the context of Riemannian geometry. For our 
purposes, the main point is that it is expected that the mean curvature blows up as the 
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singularity is approached. By using inverse mean curvature flow, the lapse gets smaller 
the nearer we approach the singularity. Inverse mean curvature flow is a parabolic 
equation, as can be seen in our case by the presence of second spatial derivatives on the 
right hand side of equation (JHj). In fact, this is the only place where second derivatives 
occur, and the entire evolution is of a mixed hyperbolic-parabolic nature. [10] The other 
nice property of the system has to do with the use of scale invariant variables: thus for 
example, as the singularity is approached, one would expect both the mean curvature 
3H and the shear a a/ 3 to blow up, but the scale invariant quantity S a/ g = o- a p/H to 
remain finite. 

Nonetheless, a mixed hyperbolic-parabolic system is somewhat unusual from the 
point of view of numerical relativity. Essentially it comes about because while we expect 
gravity to be described by a hyperbolic system, we always have a choice of gauge, and in 
this particular case we have chosen one that is parabolic. This is somewhat analogous to 
a more usual situation in numerical relativity where one chooses a gauge like maximal 
slicing which is an elliptic condition, and the vacuum Einstein equation then becomes 
a mixed hyperbolic-elliptic system. Lars Andersson has suggested that one could study 
singularities by simply using the system of ref. [8] but with constant mean curvature 
time slices as the gauge instead of inverse mean curvature flow, [llj To implement this 
suggestion, we simply choose the time coordinate t so that 3H = e~*. (Note that the 
singularity is approached as t — ► — oo). This gives rise to the following elliptic equation 



where the scale invariant lapse Af is defined by M = NH. Since H is now constant 
on the surfaces of constant time, it follows that r a = 0. Therefore, equation f[20l 
simply becomes q = (l/3)E Q/3 E Q/ g. Thus, the only quantities that we need to evolve 
are E a \ A a , N a/3 and S Q/ g. These quantities are subject to the following evolution 
equations. 



d a d a Af + 2A a d a M + (3 + £' 



T> a0 )M = 3 



(23) 




(24) 
(25) 



(26) 



+ Af [-3£ Q/3 - d <a A p> - 2iV <a ^iV /3>7 + N\N <aP> 
+ e, 5(a (d' l N^-2Am^)} 



(27) 



and to the vanishing of the following constraints. 




(28) 
(29) 
(30) 



C G 



= 1 + 



\d a A a - A a A a - \N aP N' 




(31) 
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To distinguish these two versions of the scale invariant tetrad system, we will call the 
first one the inverse mean curvature system and the second one the constant mean 
curvature (or just CMC) system. 

In order to evolve either one of these scale invariant systems, one must find initial 
data that satisfies the constraint equations. This is done essentially by using the York 
method to find initial data for the usual spatial metric and extrinsic curvature variables 
and then translating that initial data into data for the scale invariant tetrad variables. In 
particular, choose spatial topology T 3 , constant mean curvature, and choose the spatial 
metric to be of the form h a b = ip 4 S a b where 6 a b is a flat metric. Then, the momentum 
constraint becomes a differential equation in the flat space coordinates for the trace- 
free part of the extrinsic curvature (multiplied by a power of the conformal factor), 
which on T 3 with a conformally flat metric simply becomes an algebraic equation for 
the Fourier components of that quantity. The Hamiltonian constraint then becomes an 
elliptic equation for the conformal factor if). A simple family of initial data of this kind 
is the following: H = constant, N a ? = 0, Ej = (ip~ 2 /H)Sj , A a = -2ip~ l d a ifj and 

= ( — V ,_6 /-f^)diag(Si, S 2 , S3), with if) satisfying the equation 



with the a; and 6j constants. This set of conditions is enough to satisfy the constraints of 
the CMC system. If one wants to satisfy the constraints of the inverse mean curvature 
system, one must add the condition r a = and one must solve equation (TSUI) for q. 

We now consider the form of the BKL conjecture for these scale invariant systems. 
First note that all spatial derivatives occur in the form d a = E a l di. Thus one would 
expect spatial derivatives to be negligible as the singularity is approached provided 
that E a l becomes negligible and provided that any derivatives with respect to spatial 
coordinates at least do not become so large to overwhelm the smallness of E a \ From 
the form of the evolution equations, there are reasons to expect this to occur, [12j and 
these reasons also lead one to expect that A a (and in the inverse mean curvature case 
also r a ) become negligibly small. However, the ultimate test for all this is numerical: 
evolve initial data towards the singularity and see whether these quantities do indeed 
become negligible. 

Note that as in the Gowdy case spatial derivatives becoming negligible does not 
mean that the spacetime is becoming homogeneous. Instead, the dependence of the 
metric components on space can become quite large (especially at spikes); however 
since in the field equations these spatial derivatives are multiplied by a small quantity 
(e~ 2r in the Gowdy case and E a l in the scale invariant tetrad case) this means that the 
spatial derivatives make a negligible contribution to the field equations. Nonetheless, 




(32) 



and the given by 



Si = a 2 cos y + a 3 cos z + b 2 + b 3 
£ 2 = a-i cos x — a 3 cos z + 61 — b 3 
£3 = — a\ cos x — a 2 cos y — b\ — b 2 



(33) 
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with spatial derivatives becoming negligible in the field equations, this means that 
the dynamics at each spatial point is that of a homogeneous spacetime (albeit a 
different homogeneous spacetime for each spatial point). Since there are many different 
homogeneous spacetimes (classified according to the Bianchi classification) there remains 
the question of which homogeneous spacetimes describe the dynamics of the general 
singularity. Another part of the BKL conjecture is that it is the dynamics of the Bianchi 
type IX spacetimes, also called Mixmaster dynamics, that describes the approach to 
the singularity. This dynamics consists of "epochs" where the components of E Q( g are 
approximately constant, punctuated by short "bounces" where they change rapidly. The 
components of N a /3 are negligible during the epochs and exhibit rapid growth and equally 
rapid decay during the bounces. All this should also be checked numerically. However 
if the spatial derivative terms, as well as A a (and r a in the inverse mean curvature case) 
are negligible, then one can replace the evolution equations and constraint equations 
with a set of truncated equations where these terms are absent. An analysis of these 
truncated equations [12] (which due to the neglect of the spatial derivatives are ODEs) 
shows that this sort of local Mixmaster dynamics is exactly what is to be expected. Thus 
the main thing to be checked from the numerics is whether the terms that are expected 
to become negligible do in fact become negligible as the singularity is approached. 

We now turn to numerical methods and presentation of numerical results. The 
results for the inverse mean curvature system have been presented in. [12 so here we 
will concentrate on the CMC system. The spatial topology is chosen to be T 3 with 
coordinates (x, y, z) satisfying < x, y, z < 2n and with and 2n identified. This 
is implemented numerically by using n + 2 grid points in each spatial direction with 
spacing dx = 2Ti/n between adjacent grid points. The variables on gridpoints 2 to n+ 1 
are evolved, while those at gridpoints 1 and n + 2 are set using the periodic boundary 
conditions. Derivatives with respect to spatial coordinates are evaluated using centered 
differences, while time derivatives are implemented using the ICN method. Equation 
( |23|) for the scale invariant lapse is solved at each time step using a multigrid method. 
The initial data is of the form in equation (1331) with equation (1321) for the conformal 
factor solved by relaxation. 

The simulation was done on a SunBlade 2000 with n = 32 and with H = 1/3, = 
(0.2,0.1,0.02) and hi = (1.8,-0.15,0). Figure [3] shows for each time the logarithms of 
the maximum values (over all spatial points and all indicies) of \E a l \ and \A a \. Note 
that these quantities become negligibly small and so the BKL conjecture is supported. 
Since spatial derivatives are negligibly small, the dynamics is local and can be studied 
by simply looking at the behavior with time of the scale invariant quantities at a single 
spatial point. Figures H] and [5] show this sort of behavior for the diagonal components 
of E Q/ g and N a p. Note that sufficiently near the singularity, the behavior is a series 
of epochs where the components of E Qj g are approximately constant, punctuated by 
bounces where the components of E Q/3 change rapidly. Also note that N a/3 is negligible 
except at the bounces. The tensor E a) g has three eigenvalues. However, E^ is traceless, 
and during each epoch satisfies E a/3 E Q/3 = 6. The eigenvalues of can therefore be 
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Figure 5. Plot of the components of N a p vs t 




Figure 6. Plot of u vs t 
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determined by a single quantity u defined by 




(34) 



(There is a unique u > 1 provided that the quantity on the left hand side of the equation 
is between —6 and 6). It is a consequence of Mixmaster dynamics that u is constant in 
each epoch and that from one epoch to the next u changes by the map u — > u — 1 if 
u > 2 and u — > 1/ (it — 1) ifl < w < 2. figure [6] shows the value of w at the same spatial 
point. For times sufficiently near the singularity the sequence of u values satisfies this 
map. 

4. open problems 

The main limitation of the results just presented is the low resolution of the simulation 
(though simulations of this sort of system but assuming symmetry have been done with 
high resolution 1 13]). This is especially true because the general singularity has spikes 
analogous to those of the Gowdy spacetime. In the system of scale invariant tetrad 
variables, these occur because of the nature of the bounces: a bounce occurs because in 
any epoch one component of N a p grows exponentially with time until it becomes large 
enough to cause a sudden change in the components of S a/3 . This occurs whether that 
component of N a p is positive or negative; but if it is zero then it remains zero and there 
is no bounce. Thus in general spikes occur on surfaces of co-dimension 1, so what were 
isolated points in the Gowdy case are isolated surfaces in the general case. To properly 
treat the spikes in the general case will require a singularity simulating code that is 
parallel and uses adaptive mesh refinement. Development of such a code is presently 
work in progress. 

Another limitation of the current results is that they are only for the vacuum case. 
Another part of the BKL conjecture is that as the singularity is approached, matter 
makes a negligible contribution to the gravitational field equations unless it is a stiff 
fluid (equation of state P = p). In the stiff fluid case, BKL conjecture that there is a 
last bounce after which the components of H a p remain essentially constant all the way 
to the singularity. This stiff fluid part of the BKL conjecture is supported by a theorem 
of Rendall and Andersson[14] as well as by numerical simulations. [15] What is needed is 
some numerical simulations to cover the non-stiff fluid case. In (8] scale invariant tetrad 
equations are given for fluids with equation of state P = kp for constant k. What is 
needed is a numerical simulation of these equations to see whether for k < 1 the fluid 
has a negligible influence on the gravitational degrees of freedom as the singualarity is 
approached. Since fluids with k < 1 can form shocks, such a numerical simulation would 
need a high resolution shock capturing method. 

Another limitation of the current results is that they are for spacetimes with 
compact spatial topology and therefore do not directly address the physically relevant 
case of gravitational collapse in an asymptotically flat spacetime. This is not as bad as 
it seems: since the simulations show that the formation of singularities in spacetimes 
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with T 3 spatial topology is essentially a local process, this means that one can expect 
at least parts of the singularities produced in asymptotically flat spacetimes to have the 
same features as the singuarities in these simulations. However, that still leaves open the 
possibility that for gravitational collapse in asymptotically flat spacetimes other parts 
of the singularity have very different properties. In particularly Poisson and Israel have 
conjectured that the part of the singularity near the horizon is not of the BKL type 
but is instead null rather than spacelike and with weaker tidal forces. [16] Numerical 
simulations in the spherically symmetric case support this conjecture. |17[ fT8] What is 
needed now is numerical simulations that don't assume spherical symmetry. 
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